step2Model2ComputingEff_2012_Inc = function(sims){
  ATEm20=rep(NA, sims)
  AMEm20=rep(NA, sims)
  ADEm20=rep(NA, sims)
  ATEm15=rep(NA, sims)
  AMEm15=rep(NA, sims)
  ADEm15=rep(NA, sims)
  ATEm10=rep(NA, sims)
  AMEm10=rep(NA, sims)
  ADEm10=rep(NA, sims)
  ATEm5=rep(NA, sims)
  AMEm5=rep(NA, sims)
  ADEm5=rep(NA, sims)
  ATE0=rep(NA, sims)
  AME0=rep(NA, sims)
  ADE0=rep(NA, sims)
  ATE5=rep(NA, sims)
  AME5=rep(NA, sims)
  ADE5=rep(NA, sims)
  ATE10=rep(NA, sims)
  AME10=rep(NA, sims)
  ADE10=rep(NA, sims)
  ATE15=rep(NA, sims)
  AME15=rep(NA, sims)
  ADE15=rep(NA, sims)
  ATE20=rep(NA, sims)
  AME20=rep(NA, sims)
  ADE20=rep(NA, sims)
  ATE25=rep(NA, sims)
  AME25=rep(NA, sims)
  ADE25=rep(NA, sims)
  
  for (i in 1:sims){
    ## TOTAL EFFECT ##
    newxm25 = cbind(t(holder[1:5,,(i-1)*11+1]),
                     df$chgUnemp1211, df$avgUnemp12, rep(-2.5, nrow(df)), df$inc12, df$avgGas, df$totForcAugOct, 
                     df$younger30, df$older65, df$female, df$black, df$latino, df$noHighSchool, df$HighSchool, df$someCollege, df$fourCollege,  
                     df$income, df$unemployed, df$fullTime, df$partTime, df$ownHome,
                     df$Very, df$Few, df$ideol, df$Dem, df$Rep, df$Ind)
    newxm20 = cbind(t(holder[1:5,,(i-1)*11+2]),
                     df$chgUnemp1211, df$avgUnemp12, rep(-2.0, nrow(df)), df$inc12, df$avgGas, df$totForcAugOct, 
                     df$younger30, df$older65, df$female, df$black, df$latino, df$noHighSchool, df$HighSchool, df$someCollege, df$fourCollege,  
                     df$income, df$unemployed, df$fullTime, df$partTime, df$ownHome,
                     df$Very, df$Few, df$ideol, df$Dem, df$Rep, df$Ind)
    newxm15 = cbind(t(holder[1:5,,(i-1)*11+3]),
                     df$chgUnemp1211, df$avgUnemp12, rep(-1.5, nrow(df)), df$inc12, df$avgGas, df$totForcAugOct, 
                     df$younger30, df$older65, df$female, df$black, df$latino, df$noHighSchool, df$HighSchool, df$someCollege, df$fourCollege,  
                     df$income, df$unemployed, df$fullTime, df$partTime, df$ownHome,
                     df$Very, df$Few, df$ideol, df$Dem, df$Rep, df$Ind)
    newxm10 = cbind(t(holder[1:5,,(i-1)*11+4]),
                     df$chgUnemp1211, df$avgUnemp12, rep(-1.0, nrow(df)), df$inc12, df$avgGas, df$totForcAugOct, 
                     df$younger30, df$older65, df$female, df$black, df$latino, df$noHighSchool, df$HighSchool, df$someCollege, df$fourCollege,  
                     df$income, df$unemployed, df$fullTime, df$partTime, df$ownHome,
                     df$Very, df$Few, df$ideol, df$Dem, df$Rep, df$Ind)
    newxm5= cbind(t(holder[1:5,,(i-1)*11+5]),
                   df$chgUnemp1211, df$avgUnemp12, rep(-.5, nrow(df)), df$inc12, df$avgGas, df$totForcAugOct, 
                   df$younger30, df$older65, df$female, df$black, df$latino, df$noHighSchool, df$HighSchool, df$someCollege, df$fourCollege,  
                   df$income, df$unemployed, df$fullTime, df$partTime, df$ownHome,
                   df$Very, df$Few, df$ideol, df$Dem, df$Rep, df$Ind)
    newx0 = cbind(t(holder[1:5,,(i-1)*11+6]),
                   df$chgUnemp1211, df$avgUnemp12, rep(0, nrow(df)), df$inc12, df$avgGas, df$totForcAugOct, 
                   df$younger30, df$older65, df$female, df$black, df$latino, df$noHighSchool, df$HighSchool, df$someCollege, df$fourCollege,  
                   df$income, df$unemployed, df$fullTime, df$partTime, df$ownHome,
                   df$Very, df$Few, df$ideol, df$Dem, df$Rep, df$Ind)
    newx5 = cbind(t(holder[1:5,,(i-1)*11+7]),
                   df$chgUnemp1211, df$avgUnemp12, rep(.5, nrow(df)), df$inc12, df$avgGas, df$totForcAugOct, 
                   df$younger30, df$older65, df$female, df$black, df$latino, df$noHighSchool, df$HighSchool, df$someCollege, df$fourCollege,  
                   df$income, df$unemployed, df$fullTime, df$partTime, df$ownHome,
                   df$Very, df$Few, df$ideol, df$Dem, df$Rep, df$Ind)
    newx10 = cbind(t(holder[1:5,,(i-1)*11+8]),
                    df$chgUnemp1211, df$avgUnemp12, rep(1.0, nrow(df)), df$inc12, df$avgGas, df$totForcAugOct, 
                    df$younger30, df$older65, df$female, df$black, df$latino, df$noHighSchool, df$HighSchool, df$someCollege, df$fourCollege,  
                    df$income, df$unemployed, df$fullTime, df$partTime, df$ownHome,
                    df$Very, df$Few, df$ideol, df$Dem, df$Rep, df$Ind)
    newx15 = cbind(t(holder[1:5,,(i-1)*11+9]),
                    df$chgUnemp1211, df$avgUnemp12, rep(1.5, nrow(df)), df$inc12, df$avgGas, df$totForcAugOct, 
                    df$younger30, df$older65, df$female, df$black, df$latino, df$noHighSchool, df$HighSchool, df$someCollege, df$fourCollege,  
                    df$income, df$unemployed, df$fullTime, df$partTime, df$ownHome,
                    df$Very, df$Few, df$ideol, df$Dem, df$Rep, df$Ind)
    newx20 = cbind(t(holder[1:5,,(i-1)*11+10]),
                    df$chgUnemp1211, df$avgUnemp12, rep(2.0, nrow(df)), df$inc12, df$avgGas, df$totForcAugOct, 
                    df$younger30, df$older65, df$female, df$black, df$latino, df$noHighSchool, df$HighSchool, df$someCollege, df$fourCollege,  
                    df$income, df$unemployed, df$fullTime, df$partTime, df$ownHome,
                    df$Very, df$Few, df$ideol, df$Dem, df$Rep, df$Ind)
    newx25 = cbind(t(holder[1:5,,(i-1)*11+11]),
                    df$chgUnemp1211, df$avgUnemp12, rep(2.5, nrow(df)), df$inc12, df$avgGas, df$totForcAugOct, 
                    df$younger30, df$older65, df$female, df$black, df$latino, df$noHighSchool, df$HighSchool, df$someCollege, df$fourCollege,  
                    df$income, df$unemployed, df$fullTime, df$partTime, df$ownHome,
                    df$Very, df$Few, df$ideol, df$Dem, df$Rep, df$Ind)
    
    ## MEDIATION EFFECT ##
    newxm25M = cbind(t(holder[1:5,,(i-1)*11+1]),
                      df$chgUnemp1211, df$avgUnemp12, rep(-1.5, nrow(df)), df$inc12, df$avgGas, df$totForcAugOct, 
                      df$younger30, df$older65, df$female, df$black, df$latino, df$noHighSchool, df$HighSchool, df$someCollege, df$fourCollege,  
                      df$income, df$unemployed, df$fullTime, df$partTime, df$ownHome,
                      df$Very, df$Few, df$ideol, df$Dem, df$Rep, df$Ind)
    newxm20M = cbind(t(holder[1:5,,(i-1)*11+2]),
                      df$chgUnemp1211, df$avgUnemp12, rep(-1.0, nrow(df)), df$inc12, df$avgGas, df$totForcAugOct, 
                      df$younger30, df$older65, df$female, df$black, df$latino, df$noHighSchool, df$HighSchool, df$someCollege, df$fourCollege,  
                      df$income, df$unemployed, df$fullTime, df$partTime, df$ownHome,
                      df$Very, df$Few, df$ideol, df$Dem, df$Rep, df$Ind)
    newxm15M = cbind(t(holder[1:5,,(i-1)*11+3]), 
                      df$chgUnemp1211, df$avgUnemp12, rep(-.5, nrow(df)), df$inc12, df$avgGas, df$totForcAugOct, 
                      df$younger30, df$older65, df$female, df$black, df$latino, df$noHighSchool, df$HighSchool, df$someCollege, df$fourCollege,  
                      df$income, df$unemployed, df$fullTime, df$partTime, df$ownHome,
                      df$Very, df$Few, df$ideol, df$Dem, df$Rep, df$Ind)
    newxm10M = cbind(t(holder[1:5,,(i-1)*11+4]), 
                      df$chgUnemp1211, df$avgUnemp12, rep(0, nrow(df)), df$inc12, df$avgGas, df$totForcAugOct, 
                      df$younger30, df$older65, df$female, df$black, df$latino, df$noHighSchool, df$HighSchool, df$someCollege, df$fourCollege,  
                      df$income, df$unemployed, df$fullTime, df$partTime, df$ownHome,
                      df$Very, df$Few, df$ideol, df$Dem, df$Rep, df$Ind)
    newxm5M = cbind(t(holder[1:5,,(i-1)*11+5]), 
                     df$chgUnemp1211, df$avgUnemp12, rep(.5, nrow(df)), df$inc12, df$avgGas, df$totForcAugOct, 
                     df$younger30, df$older65, df$female, df$black, df$latino, df$noHighSchool, df$HighSchool, df$someCollege, df$fourCollege,  
                     df$income, df$unemployed, df$fullTime, df$partTime, df$ownHome,
                     df$Very, df$Few, df$ideol, df$Dem, df$Rep, df$Ind)
    newx0M = cbind(t(holder[1:5,,(i-1)*11+6]), 
                    df$chgUnemp1211, df$avgUnemp12, rep(1.0, nrow(df)), df$inc12, df$avgGas, df$totForcAugOct, 
                    df$younger30, df$older65, df$female, df$black, df$latino, df$noHighSchool, df$HighSchool, df$someCollege, df$fourCollege,  
                    df$income, df$unemployed, df$fullTime, df$partTime, df$ownHome,
                    df$Very, df$Few, df$ideol, df$Dem, df$Rep, df$Ind)
    newx5M = cbind(t(holder[1:5,,(i-1)*11+7]), 
                    df$chgUnemp1211, df$avgUnemp12, rep(1.5, nrow(df)), df$inc12, df$avgGas, df$totForcAugOct, 
                    df$younger30, df$older65, df$female, df$black, df$latino, df$noHighSchool, df$HighSchool, df$someCollege, df$fourCollege,  
                    df$income, df$unemployed, df$fullTime, df$partTime, df$ownHome,
                    df$Very, df$Few, df$ideol, df$Dem, df$Rep, df$Ind)
    newx10M = cbind(t(holder[1:5,,(i-1)*11+8]), 
                     df$chgUnemp1211, df$avgUnemp12, rep(2.0, nrow(df)), df$inc12, df$avgGas, df$totForcAugOct, 
                     df$younger30, df$older65, df$female, df$black, df$latino, df$noHighSchool, df$HighSchool, df$someCollege, df$fourCollege,  
                     df$income, df$unemployed, df$fullTime, df$partTime, df$ownHome,
                     df$Very, df$Few, df$ideol, df$Dem, df$Rep, df$Ind)
    newx15M = cbind(t(holder[1:5,,(i-1)*11+9]), 
                     df$chgUnemp1211, df$avgUnemp12, rep(2.5, nrow(df)), df$inc12, df$avgGas, df$totForcAugOct, 
                     df$younger30, df$older65, df$female, df$black, df$latino, df$noHighSchool, df$HighSchool, df$someCollege, df$fourCollege,  
                     df$income, df$unemployed, df$fullTime, df$partTime, df$ownHome,
                     df$Very, df$Few, df$ideol, df$Dem, df$Rep, df$Ind)
    
    ## DIRECT EFFECT ##
    newxm15D = cbind(t(holder[1:5,,(i-1)*11+1]),
                      df$chgUnemp1211, df$avgUnemp12, rep(-1.5, nrow(df)), df$inc12, df$avgGas, df$totForcAugOct, 
                      df$younger30, df$older65, df$female, df$black, df$latino, df$noHighSchool, df$HighSchool, df$someCollege, df$fourCollege,  
                      df$income, df$unemployed, df$fullTime, df$partTime, df$ownHome,
                      df$Very, df$Few, df$ideol, df$Dem, df$Rep, df$Ind)
    newxm10D = cbind(t(holder[1:5,,(i-1)*11+2]),
                      df$chgUnemp1211, df$avgUnemp12, rep(-1.0, nrow(df)), df$inc12, df$avgGas, df$totForcAugOct, 
                      df$younger30, df$older65, df$female, df$black, df$latino, df$noHighSchool, df$HighSchool, df$someCollege, df$fourCollege,  
                      df$income, df$unemployed, df$fullTime, df$partTime, df$ownHome,
                      df$Very, df$Few, df$ideol, df$Dem, df$Rep, df$Ind)
    newxm5D = cbind(t(holder[1:5,,(i-1)*11+3]), 
                     df$chgUnemp1211, df$avgUnemp12, rep(-.5, nrow(df)), df$inc12, df$avgGas, df$totForcAugOct, 
                     df$younger30, df$older65, df$female, df$black, df$latino, df$noHighSchool, df$HighSchool, df$someCollege, df$fourCollege,  
                     df$income, df$unemployed, df$fullTime, df$partTime, df$ownHome,
                     df$Very, df$Few, df$ideol, df$Dem, df$Rep, df$Ind)
    newx0D = cbind(t(holder[1:5,,(i-1)*11+4]), 
                    df$chgUnemp1211,df$avgUnemp12, rep(0, nrow(df)), df$inc12, df$avgGas, df$totForcAugOct, 
                    df$younger30, df$older65, df$female, df$black, df$latino, df$noHighSchool, df$HighSchool, df$someCollege, df$fourCollege,  
                    df$income, df$unemployed, df$fullTime, df$partTime, df$ownHome,
                    df$Very, df$Few, df$ideol, df$Dem, df$Rep, df$Ind)
    newx5D = cbind(t(holder[1:5,,(i-1)*11+5]), 
                    df$chgUnemp1211, df$avgUnemp12, rep(.5, nrow(df)), df$inc12, df$avgGas, df$totForcAugOct, 
                    df$younger30, df$older65, df$female, df$black, df$latino, df$noHighSchool, df$HighSchool, df$someCollege, df$fourCollege,  
                    df$income, df$unemployed, df$fullTime, df$partTime, df$ownHome,
                    df$Very, df$Few, df$ideol, df$Dem, df$Rep, df$Ind)
    newx10D = cbind(t(holder[1:5,,(i-1)*11+6]), 
                     df$chgUnemp1211, df$avgUnemp12, rep(1.0, nrow(df)), df$inc12, df$avgGas, df$totForcAugOct, 
                     df$younger30, df$older65, df$female, df$black, df$latino, df$noHighSchool, df$HighSchool, df$someCollege, df$fourCollege,  
                     df$income, df$unemployed, df$fullTime, df$partTime, df$ownHome,
                     df$Very, df$Few, df$ideol, df$Dem, df$Rep, df$Ind)
    newx15D = cbind(t(holder[1:5,,(i-1)*11+7]), 
                     df$chgUnemp1211, df$avgUnemp12, rep(1.5, nrow(df)), df$inc12, df$avgGas, df$totForcAugOct, 
                     df$younger30, df$older65, df$female, df$black, df$latino, df$noHighSchool, df$HighSchool, df$someCollege, df$fourCollege,  
                     df$income, df$unemployed, df$fullTime, df$partTime, df$ownHome,
                     df$Very, df$Few, df$ideol, df$Dem, df$Rep, df$Ind)
    newx20D = cbind(t(holder[1:5,,(i-1)*11+8]), 
                     df$chgUnemp1211, df$avgUnemp12, rep(2.0, nrow(df)), df$inc12, df$avgGas, df$totForcAugOct, 
                     df$younger30, df$older65, df$female, df$black, df$latino, df$noHighSchool, df$HighSchool, df$someCollege, df$fourCollege,  
                     df$income, df$unemployed, df$fullTime, df$partTime, df$ownHome,
                     df$Very, df$Few, df$ideol, df$Dem, df$Rep, df$Ind)
    newx25D = cbind(t(holder[1:5,,(i-1)*11+9]), 
                     df$chgUnemp1211, df$avgUnemp12, rep(2.5, nrow(df)), df$inc12, df$avgGas, df$totForcAugOct, 
                     df$younger30, df$older65, df$female, df$black, df$latino, df$noHighSchool, df$HighSchool, df$someCollege, df$fourCollege,  
                     df$income, df$unemployed, df$fullTime, df$partTime, df$ownHome,
                     df$Very, df$Few, df$ideol, df$Dem, df$Rep, df$Ind)
    
    
    
    coeffics = matrix(as.numeric(maincoefs[i,]), ncol = 1)
    evalRandomState = matrix(NA, nrow=length(groupcoefsID), ncol=2) # Create an empty matrix that will store random effects
    evalRandomState[,1] = t(groupcoefs[i, ]) # Each row of groupcoefs has intercepts
    evalRandomState[,2]=groupcoefsID[1:length(groupcoefsID)] # Col 2 is for stateID
    colnames(evalRandomState) = c("Intercept", "stateID")
    newf=merge(grouponly, evalRandomState, by = "stateID", all.x = TRUE)
    newf = as.data.frame(cbind(newf$Intercept, newf$IDnum))
    names(newf) = c("Intercept", "IDnum")
    newf = newf[order(newf$IDnum),]
    
    etam25 = newxm25%*%coeffics + newf[,1] # newf[,1]: random effects
    etam20 = newxm20%*%coeffics + newf[,1]
    etam15 = newxm15%*%coeffics + newf[,1]
    etam10 = newxm10%*%coeffics + newf[,1]
    etam5 = newxm5%*%coeffics + newf[,1]
    eta0 = newx0%*%coeffics + newf[,1]
    eta5 = newx5%*%coeffics + newf[,1]
    eta10 = newx10%*%coeffics + newf[,1]
    eta15 = newx15%*%coeffics + newf[,1]
    eta20 = newx20%*%coeffics + newf[,1]
    eta25 = newx25%*%coeffics + newf[,1]
    
    etaAMEm25 = newxm25M%*%coeffics + newf[,1]
    etaAMEm20 = newxm20M%*%coeffics + newf[,1]
    etaAMEm15 = newxm15M%*%coeffics + newf[,1]
    etaAMEm10 = newxm10M%*%coeffics + newf[,1]
    etaAMEm5 = newxm5M%*%coeffics + newf[,1]
    etaAME0 = newx0M%*%coeffics + newf[,1]
    etaAME5 = newx5M%*%coeffics + newf[,1]
    etaAME10 = newx10M%*%coeffics + newf[,1]
    etaAME15 = newx15M%*%coeffics + newf[,1]
    
    etaADEm15 = newxm15D%*%coeffics + newf[,1]
    etaADEm10 = newxm10D%*%coeffics + newf[,1]
    etaADEm5 = newxm5D%*%coeffics + newf[,1]
    etaADE0 = newx0D%*%coeffics + newf[,1]
    etaADE5 = newx5D%*%coeffics + newf[,1]
    etaADE10 = newx10D%*%coeffics + newf[,1]
    etaADE15 = newx15D%*%coeffics + newf[,1]
    etaADE20 = newx20D%*%coeffics + newf[,1]
    etaADE25 = newx25D%*%coeffics + newf[,1]
    
    
    ATEm15[i] = mean(exp(etam15)/(1 + exp(etam15)), na.rm = T) - mean(exp(etam25)/(1 + exp(etam25)), na.rm = T)
    AMEm15[i] = mean(exp(etam15)/(1 + exp(etam15)), na.rm = T) - mean(exp(etaAMEm25)/(1 + exp(etaAMEm25)), na.rm = T)
    ADEm15[i] = mean(exp(etaADEm15)/(1 + exp(etaADEm15)), na.rm = T) - mean(exp(etam25)/(1 + exp(etam25)), na.rm = T)
    
    ATEm10[i] = mean(exp(etam10)/(1 + exp(etam10)), na.rm = T) - mean(exp(etam20)/(1 + exp(etam20)), na.rm = T)
    AMEm10[i] = mean(exp(etam10)/(1 + exp(etam10)), na.rm = T) - mean(exp(etaAMEm20)/(1 + exp(etaAMEm20)), na.rm = T)
    ADEm10[i] = mean(exp(etaADEm10)/(1 + exp(etaADEm10)), na.rm = T) - mean(exp(etam20)/(1 + exp(etam20)), na.rm = T)
    
    ATEm5[i] = mean(exp(etam5)/(1 + exp(etam5)), na.rm = T) - mean(exp(etam15)/(1 + exp(etam15)), na.rm = T)
    AMEm5[i] = mean(exp(etam5)/(1 + exp(etam5)), na.rm = T) - mean(exp(etaAMEm15)/(1 + exp(etaAMEm15)), na.rm = T)
    ADEm5[i] = mean(exp(etaADEm5)/(1 + exp(etaADEm5)), na.rm = T) - mean(exp(etam15)/(1 + exp(etam15)), na.rm = T)
    
    ATE0[i] = mean(exp(eta0)/(1 + exp(eta0)), na.rm = T) - mean(exp(etam10)/(1 + exp(etam10)), na.rm = T)
    AME0[i] = mean(exp(eta0)/(1 + exp(eta0)), na.rm = T) - mean(exp(etaAMEm10)/(1 + exp(etaAMEm10)), na.rm = T)
    ADE0[i] = mean(exp(etaADE0)/(1 + exp(etaADE0)), na.rm = T) - mean(exp(etam10)/(1 + exp(etam10)), na.rm = T)
    
    ATE5[i] = mean(exp(eta5)/(1 + exp(eta5)), na.rm = T) - mean(exp(etam5)/(1 + exp(etam5)), na.rm = T)
    AME5[i] = mean(exp(eta5)/(1 + exp(eta5)), na.rm = T) - mean(exp(etaAMEm5)/(1 + exp(etaAMEm5)), na.rm = T)
    ADE5[i] = mean(exp(etaADE5)/(1 + exp(etaADE5)), na.rm = T) - mean(exp(etam5)/(1 + exp(etam5)), na.rm = T)
    
    ATE10[i] = mean(exp(eta10)/(1 + exp(eta10)), na.rm = T) - mean(exp(eta0)/(1 + exp(eta0)), na.rm = T)
    AME10[i] = mean(exp(eta10)/(1 + exp(eta10)), na.rm = T) - mean(exp(etaAME0)/(1 + exp(etaAME0)), na.rm = T)
    ADE10[i] = mean(exp(etaADE10)/(1 + exp(etaADE10)), na.rm = T) - mean(exp(eta0)/(1 + exp(eta0)), na.rm = T)
    
    ATE15[i] = mean(exp(eta15)/(1 + exp(eta15)), na.rm = T) - mean(exp(eta5)/(1 + exp(eta5)), na.rm = T)
    AME15[i] = mean(exp(eta15)/(1 + exp(eta15)), na.rm = T) - mean(exp(etaAME5)/(1 + exp(etaAME5)), na.rm = T)
    ADE15[i] = mean(exp(etaADE15)/(1 + exp(etaADE15)), na.rm = T) - mean(exp(eta5)/(1 + exp(eta5)), na.rm = T)
    
    ATE20[i] = mean(exp(eta20)/(1 + exp(eta20)), na.rm = T) - mean(exp(eta10)/(1 + exp(eta10)), na.rm = T)
    AME20[i] = mean(exp(eta20)/(1 + exp(eta20)), na.rm = T) - mean(exp(etaAME10)/(1 + exp(etaAME10)), na.rm = T)
    ADE20[i] = mean(exp(etaADE20)/(1 + exp(etaADE20)), na.rm = T) - mean(exp(eta10)/(1 + exp(eta10)), na.rm = T)
    
    ATE25[i] = mean(exp(eta25)/(1 + exp(eta25)), na.rm = T) - mean(exp(eta15)/(1 + exp(eta15)), na.rm = T)
    AME25[i] = mean(exp(eta25)/(1 + exp(eta25)), na.rm = T) - mean(exp(etaAME15)/(1 + exp(etaAME15)), na.rm = T)
    ADE25[i] = mean(exp(etaADE25)/(1 + exp(etaADE25)), na.rm = T) - mean(exp(eta15)/(1 + exp(eta15)), na.rm = T)
    
    if (i %% 20 == 0) print(i)
    if (i %% 20 == 0) cat("Will be done when this number becomes", sims, ". Keep cool :) \n") 
  }
  return(as.data.frame(cbind(ATEm15, AMEm15, ADEm15, 
                             ATEm10, AMEm10, ADEm10, 
                             ATEm5, AMEm5, ADEm5, 
                             ATE0, AME0, ADE0, 
                             ATE5, AME5, ADE5, 
                             ATE10, AME10, ADE10,
                             ATE15, AME15, ADE15,
                             ATE20, AME20, ADE20,
                             ATE25, AME25, ADE25)))
}